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ABSTRACT 

We develop a simple evolutionary scenario for the growth of supermassive black holes (BHs), 
assuming growth due to accretion only, to learn about the evolution of the BH mass function 
from z = 3 to and from it calculate the energy budgets of different modes of feedback. 
We tune the parameters of the model by matching the derived X-ray luminosity function 
(XLF) with the observed XLF of active galactic nuclei. We then calculate the amount of co- 
moving kinetic and bolometric feedback as a function of redshift, derive a kinetic luminosity 
function and estimate the amount of kinetic feedback and PdV work done by classical dou- 
ble Fanaroff-Riley II (FR II) radio sources. We also derive the radio luminosity function for 
FR lis from our synthesized population and set constraints on jet duty cycles. Around 1/6 of 
the jet power from FR II sources goes into PdV work done in the expanding lobes during the 
time the jet is on. Anti hierarchical growth of BHs is seen in our model due to addition of an 
amount of mass being accreted on to all BHs independent of the BH mass. The contribution 
to the total kinetic feedback by active galaxies in a low accretion, kinetically efficient mode is 
found to be the most significant at z < 1.5. FR II feedback is found to be a significant mode 
of feedback above redshifts z ~ 1.5, which has not been highlighted by previous studies. 
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1 INTRODUCTION 

As supermassive black holes (BHs) at centres of galaxies grow 
through mass accretion, they also exhibit different types of feed- 
back in radiative and kinetic forms. Understanding the energy 
budget of the various feedback components is important because 
they may couple to the surrounding environment in distinct ways. 
Most, if not all, galaxi es are believed to host supermassive BHs 
(Magorrian et al. 1998), and their roles in galaxy evolution can- 
not be neglected. Mechanical feedback from active galactic nuclei 
(AGN) has recently been invoked to be an important part of galaxy 
formation and the heating of inflowing gas, providing energy to 
suppress star f ormation (first steps have been taken by a num- 
ber of groups: Kauffmann & Haehnelt 2000; Croton ^t alj|2006l : 
Hopkins & Hernquist 2006: lBower et alj|200^) . Without including 
such an additional heating mechanism, cosmologic al models over- 
predi c t the number of faint and bright galaxies (e .g. IWhite & Reed 
Il978t IWMte&Frenklll99ll : iBenson et al]|2003ri and predict that 
the largest galaxies in the present epoch are blue and star forming 
rather than 'red and dead', contrary to observations. Additionally, 
BH feedback is know n to play an impor tant role in energy input 
in galaxy clusters (e.g. iFabian etal ] l2003h and galaxy groups (e.g. 
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iMcCarthv et"ai]|2010h . Accounting for different feedback modes 
more carefully may augment our picture of galaxy evolution. 

It can be difficult to estimate the cosmic history of the en- 
ergy budget of certain feedback modes directly, though recent cal- 
culations for mechanical feedback from AGN, relying on variou s 
assumptions, have been carried out by [Merloni & Heinzl (12008 ). 
iKording. Jester & Fenderl (2008), and Icattaneo & Best! 120091) . 
One of the more challenging modes to calculate, for example, is 
mechanical feedback of classical double Fanaroff-Riley II (FR II) 
sources. Here one must carefully treat the complicated physics of 
radio emission and particle acceleration in the hotspots. The lobes 
are radio-luminous because of the synchrotron radiation from ul- 
trarelativistic electrons accelerated at the hotspots, but this may ac- 
count for only a small fraction of the total energy output, the rest of 
which may be used to do PdV work in expanding the lobes or be 
stored as internal energy. Additionally, the duty cycle of jets is not 
well known. A further complicating fa ctor is the 'youth-redshift' 
degeneracy (Blundell & Rawlings 1999), namely that older radio 
sources do not survive above previous radio survey flux limits. 
Therefore, in this work we seek an approach to estimate the cos- 
mic volume-averaged energies from different modes of feedback 
which is suited to obtain a better handle on FR II feedback than ap- 
proaches taken in previous studies. In our approach we do not make 
assumptions about jet duty cycle or relations between the total radio 
and X-ray luminosity of a source. We also do not use the observed 
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radio luminosity function (RLF) of radio galaxies directly to make 
our estimates due to the confounding effects of jet duty cycle and 
the 'youth-redshift' degeneracy. 

In our approach, we first reconstruct the cosmic history of ac- 
cretion on to supermassive BHs using a simple volume-averaged 
evolutionary scenario/prescription for accretion. Our approach is 
quite different from previous work because we assume directly a 
physically-motivated functional form for the volume-averaged ac- 
cretion rate (M) as a function of BH mass and cosmic time, with a 
limited number of free parameters to tune that affect the evolution 
in independent ways. We do expect that most galaxies were more 
active in the past, duri ng the quasar era, 1.5 < z < 3, where mos t 
radio activity is seen djarvis & Rawiingsl2000l:lBarger et al.ll200lh . 
We use observational constraints to tune the parameters: namely, 
the local BH mass function (BHMF) which has been estimated us- 
ing well-known correlations between the BH mass, bulge luminos- 
ity, an d stellar velocity dispersion of the galaxy (e.g. lMarconi et all 
12004b and the AGN luminosity function (LF) in the 2-10 keV lu- 
minosity range of 10 41,5 -1 46 ' 5 erg s" 1 , w hich has been charac- 
terized up to redshift z = 3 dUeda et al]|2003l) . 

In our evolutionary description, we assume that from z — 3 
to BH growth has occurred due to accretion only. In general, 
cosmic BH growth can be due to either BH-BH merger events 
or accretion (either secular or triggered by galaxy mergers) and 
the dominant process as a function of redshift are not fully un- 
derstood. However, observational evidence suggests that much of 
AGN evolution is secular from z = 2 to th e present epoch 
dOrban de Xivrv et al.l201 ltlKocevski et alj|2012l) . BH growth due 
to accr etion only is also enough to explain the local BH mass 
density (Fabian & Iwasawa 1999; Elvis, Risaliti & Zamorani 2002; 
lYu & Tremaind 120021 see also lKawaguchi et al]|2004l) . Such argu- 
ments imply that most BHs have an accretion efficiency of 77 ~ 0.1. 
BH growth due to mergers is only expe cted to be domin ant for high 
mass BHs (with massive dark haloes) JCattaneoi r2002). Therefore, 
it is plausible that BH growth due to merging events is negligible 
for z < 3. With this assumption, the BHMF can easily be integrated 
backw ards in time using a continuity equation dSmall & Blandford 
1992) with our prescription for (M) to obtain the history of the evo- 
lution of BHs in the centres of galaxies. However, a fully detailed 
analysis would include the significance of mergers, which are a part 
of many theoretical models for quasar evolution for triggering gas 
infall and could potentially be responsible for triggering the quasar 
peak at z = 2-3. 

To link accretion with feedback, we use a physically-based 
model for AGN based on accretion in microquasars to explain the 
ratio of kinetic and radiative released ener gy as a function of BH 
accretion rate (discussed in, for example, Merl onT"& Heinz|[2008l) . 
A BH may accrete in three different modes, depending on its ac- 
cretion rate: a low accretion, kinetic (LK) mode, a high accretion 
radiative (HR) mode (no jets), and a high accretion kinetic (HK) 
mode (powerful radio jets). After reconstructing the cosmic evo- 
lution from our prescription and observational constraints of the 
BHMF and X-ray luminosity function (XLF), we can construct an 
RLF to determine the duty cycle of powerful jets. 

There have been a number of studies that investigate the evo- 
lutionary growth of BHs. Many trace the evolution of BHs us- 
ing the LF of AGN and integrate the continuity equation to ob- 
tain a local BHMF which is matched with the observationally de- 
rived one. Alternatively, the continuity equation has been integrated 
backwards in time to obtain an evolving BHMF. Some studies as- 
sume fixed Eddington ra tios and accretion efficiencies for popula- 
tions or sub-populations dMarconi et al.l2004 ; Raimundo & Fabian 



2009), some assume redshift or mass-dependen t Eddington ra- 
tios dShankar. Weinberg & Miralda-Escu di l2009h . some assume 
a so-called 'Fundamental Plane' relation between the intrinsic 
(unbeamed) radio luminosity of the jet core, acc retion-powered 
X-ray luminosity, and BH mass at all redshifts dMerlonil 120041 : 
iMerloni & Heinzl2008l) . and others as sume a log-normal or power - 
law distribution of Eddington ratios dCao & Lill20oH ICaolboiOh . 
From these studies, an anti hierarchical growth of BHs has been 
found (meaning that most lower mass BHs formed later in cosmic 
time than higher mass BHs), with possible reversal of the d ownsiz- 
ing at z ~ 2 and higher redshifts dMerloni & Heinz 2008). Down- 
sizing was first ob servationally found using X-ray surveys in 2003 
dUeda et alj|2003l) . 

In this paper, we take a different approach by considering a 
simple evolutionary scheme for accretion on to BHs. We use the 
model to derive the evolving BHMF from the local one using the 
continuity equation, and attempt to reconstruct the observed XLF 
to learn what sort of an evolutionary scenario may be possible and 
thereby infer a plausible reason for the anti hierarchical growth. 
The real evolutionary mechanism of a single BH is, of course, more 
complex, but we can gain insight into volume-averaged average BH 
growth as a function of cosmic time and what may be responsi- 
ble for the observed anti hierarchical growth. We then calculate 
the volume-averaged kinetic feedback from jets as a function of 
redshift. We also construct a kinetic luminosity function (KLF) as 
well as an RLF f or FR II sources with analytic mo del developed for 
FR II sources in Mocz, Fabian & Blundelll J20T 1) and estimate jet 
duty cycle as a function of redshift. 

In § [2] we describe a physically motivated model for the re- 
lation between kinetically and radiatively released power and the 
accretion rate of BHs. In §[3] we describe our simple framework 
for the evolution of BHs. In §[4] we present the results of matching 
the observed XLF with our model as well as the derived popula- 
tion properties and calculation of different types of feedback. Fi- 
nally in §[5] we discuss the implications of our results. We adopt 
standard cosmological parameters S7m = 0.3, £Ia = 0.7 and 
Ho = 70 km s _1 Mpc -1 in this work. 



2 SUPERMASSIVE BLACK HOLE ACCRETION 

The accretion rate of a BH may be described in Eddington units 
as m e r/Mc 2 / L-Edd, where M is the BH mass, 77 is the accre- 
tion efficiency and I/Edd is the Eddington luminosity, Z/Edd = 
AnGMm^c/ (jt where m p is proton mass and <tt is the Thom- 
son cross-section. The accretion efficiency, 77, depends on the in- 
ner boundary conditions of the accretion flow, and in the classi- 
cal no-torque case has a value of 0.057 for a Schwarzschild BH 
(no spin) up to a value of 0.42 for a maxim ally spinning Kerr BH 
dNovikov & Thornelll973l : IMerloni & Heinz! 120081) . The accretion 
efficiency tells us the maximal amount of potential energy per unit 
time that can be extracted from the BH. The quantity M is the rate 
of mass accretion on to the BH. 

This power output of an accreting BH may be in radiative or 
kinetic form, and the efficiencies of these processes depend on both 
77 and the nature of the accretion flow. The source has bolometric 
luminosity Lboi, from which the Eddington ratio can be defined: 
A = I/boi/iEdd- The radiative efficiency of the source is: e ra d = 
Lboi/(Mc 2 ) (so that A = e ra( jMc 2 /LEdd)- In addition, BHs may 
output kinetic power with efficiency ekin = L-^i u /(Mc 2 ). If the 
kinetic power output is primarily in the form of double jets, then 



© 2013 RAS, MNRAS OOO.rjJTJT] 



AGN evolution and feedback 3 



the single jet power is Q ; e t = Z/kin/2. The kinetic output may also 
be in the form of winds felundell & Kuncic 2007|). 

An accreting BH may grow in different modes, with different 
ratios of kinetic and radiative power outputs. Here we will use a 
simple physically motivated model for modes of BH growth (see 
IMerioni & Heinj|2008t fKording. Jester & Fenderll2008i and refer- 
ences therein). The picture comes from an analogy of supermas- 
sive BHs with their less massive accreting counterparts: micro- 
quasars. We can describe the kinetic and radiative output as a func- 
tion of the accretion rate rh. Above a critical accretion rate, rh cr it, 
a BH is radiatively efficient, with A oc rh. We ado pt a value 
of merit = 3 x 10~ 2 (as in IMerioni & Heinj|2008h , at which 
A = Acrit — rn CI it . A BH accreting above the critical rate may be in 
one of two different physical states, a high accretion, kinetic (HK) 
state, i.e., one with powerful radio jets, or a high accretion, radia- 
tive (HR) state (no jets). The fraction of sources with jets switched 
on we will characterize by a jet duty cycle /, which can be a func- 
tion of cosmic time. Below m c rit, a BH will evolve in a kinetically 
efficient, radiatively inefficient mode with A oc rh 2 . This mode is 
the low accretion, kinetic (LK) mode. 

We define the kinetic luminosity to be related to the Eddington 
luminosity by the factor Akin(fh), that is: Lkm = Aki n (?h)Z/Edd. 
In the HK mode Akin (rh) oc rh and in the HR mode Akin (rh) = 
m C rit- We use the following parametrizations for A(rh) and 
Akin (hi), so that the parametrizations are differentiable and easy 
to invert: 



A(rh) = 





-6 -5 -4 -3 -2-10 1 
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Figure 1. Kinetic and radiative released power in Eddington units as a func- 
tion of BH accretion rate, rh, assuming the parametrizations of equations 
(TJ and {2}. The horizontal dotted line is at A cr i t . At low accretion rates, 
sources accrete in a radiatively inefficient (LK) mode. At high accretion 
rates, sources may accrete in either a kinetically efficient mode (HK, i.e., 
a radio-loud source with jets) or kinetically inefficient mode (HR, i.e., a 
radio-quiet source with weak/no jets). 



Akin(m) = 



+ 1 



HK mode above rh c 



HR mode above m cr it 



(2) 

The relationship between accretion rate and released power in 
kinetic and radiative form is presented in Figure [D In this picture, 
powerful radio jets are episodic events dNipo ti. Blundell & Binnevl 
120051 : [Blundell &Fabianll2(5T 1) that may happen when a source 
is accreting above the critical rate rh cr it with some duty cycle /. 
Sources with accretion rates below the critical level are assumed to 
have a high amount of kinetic output (compared to radiative out- 
put) all the time. It may be a reasonable assumption that powerful 
jetted sources behave differently than weaker jetted ones (which 
we assume are more stable and are on all the time) as the evolu- 
tion of the comoving space density of powerful radio sources is 
known to be different from that of low lumin osity radio sources 
dSmolcic et afll2009l : iMcAlpine & Jarvisll201 ll) . Also, the low lu- 
minosity Fanaroff-Riley I (FR I) sources, the less-powerful coun- 
terparts of FR lis, are kn own to reside in denser environme nts than 
FR lis out toz ~ 0.5 dAuger, Becker & Fassnachtll2008l) . which 
means that FR Is reside in environments with higher particle den- 
sity and may have a more constant fuel supply. In addition, at the 
centres of X-ray brig ht cool core clusters an FR I source is al- 
most always found dBurnsll 19901: Isiinll2009l) . indicating that less- 
powerful jet activity may be more stable and long-term processes. 



3 THE EVOLUTION OF SUPERMASSIVE BLACK 
HOLES 

Here we develop a simple physically and observationally motivated 
model for the evolution of supermassive BHs in a volume averaged 



sense. We assume that the growth and energy output of the BHs 
are due to only mass accretion from an initial cosmic time corre- 
sponding to an initial redshift, Zi up to the present, z = 0. We 
seek to reconstruct the BHMF, M 8 ^ 1 ' 1 , where M is the BH 
mass and N(M, i) is the number of BHs with mass M per unit 
comoving volume as a function of time (or redshift). Integrating 
M ^gjj over some mass range at a fixed time gives the total 
mass of BHs within the mass range per unit comoving volume. 
Performing the integral over 8N g^i'^ gives the number density of 
BHs. Assuming that the total number of BHs is constant with time 
and th at growth happens only due to accretion, the continuity equa- 
tion ( Small c&Blandford 1992) describes the evolution of the mass 
function: 



d dN{M,t) 
dt dM 



+ 



d 
dM 



dN{M,t) 



dM 



(M(M,t)) 



= (3) 



where (M(M,t)) is the mean accretion rate for all BHs of mass 
M at cosmic time t. The right-hand side of the equation would be 
equal to a source function if one wanted to describe merger events 
or the formation of new BH seeds. 

Here, we consider a simple scheme for the evolution of BHs. 
The unitless accretion rate, rh, is proportional to M /M. As a BH 
grows larger, M increases which decreases rh due to the 1/M pro- 
portionality. M is also expected to decrease with time, at least in 
an average sense, due to the depletion of available infalling matter, 
reduction of infalling material by AGN feedback, or perhaps some 
other physical mechanism that shuts down accretion (of course, an 
individual BH may have M increase with time during some periods 
in its evolution, for example a change in spin may increase the ac- 
cretion rate). We assume that at an initial redshift Zi, corresponding 
to initial lookback time tj, all the BHs are accreting brightly, with 
(rh} > 1, and accretion rates will fall exponentially with time. In 
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our evolutionary scenario, accretion near the Eddington limit and 
super-Eddington accretion are important at 2 ~ Zi, however Zi is 
chosen large enough so that super-Eddington accretion is not at all 
prevalent by the redshift 2 = 3 and older Universe, which is where 
we have observational constraints and will be testing our model. 
What happens in the model for z > 3 is not an important consid- 
eration for our results and is just an extrapolation of our model for 
z < 3. 

We assume that at time ti a BH of initial mass Mi accretes 
matter, on average, at rate 



(M(M;)) = £■ (Mi + Mc, 



(4) 



r\ cot 

where M CO nst corresponds to a constant amount of mass being ac- 
creted by BHs of all masses in a time interval r? , C Z T . The vari- 

J ' 47rGm p 

able M const is a free parameter in our model. If M const = then 
equation © corresponds to all BHs accreting with a rate (m) = 1. 
The presence of the Af const term allows low-mass and high-mass 
BHs to potentially have different accretion histories in the evolu- 
tion, which is motivated by observations that suggest anti hierar- 
chical growth of BHs. 

Next we assume that the amount of mass available to be ac- 
creted by a BH may decline with time and therefore add exponen- 
tial decay factors to the time evolution. That is, 

(M(Mi,t)) = - 47rGmp (Me^T + M const e^T) (5) 

71 CO"T 

where ri and T2 are free parameters (note that t here is the lookback 
time; t — corresponds to z = 0). 

We wish to find (M(M, i)), the average accretion rate as a 
function of lookback time and mass, from this scheme. We assume 
that all BHs at the present epoch (and at 2 < zi) of a given mass 
M had similar long-time-scale accretion histories so that they all 
originated from BHs of similar initial mass Mi at ti. That is, for a 
given M and t we solve 



(6) 



M = M; + J (M(Mi,t')> dt' 
for Mi = M ; (M, t), which yields 



Mi(M,t) 



.< 1 4.7rGm ri j — , 

1 — E-n e T i - 1 



(7) 



Thus the average accretion rate as a function of time and mass, 
needed for the continuity equation, is found: 



<M(M,t)) = (M(Mi(M,t),t)) 



(8) 



where Mi(M, t) comes from equation Q and is substituted into 
equation {3). The accretion efficiency 77 is a free parameter in the 
model. 

Knowing (M(M, t)) and a boundary condition allows us to 
numerically solve for the evolution of using the continuity 

equation (equation ([3J). Equation {8} tells us (M(M, t)). As for the 
bo undary condition, we use the 10 6 -10 10 M local BHMF derived 
bv lMarconi et alj d2004) . We set zi = 6, but will solve the BHMF 
from 2 = up to 3, where the XLF is constrained dUeda et ail 
120031) . We derive the evolving BHMF for masses 10 6 -10 10 M . 
The local BHMF is extrapolated for masses beyond 10 10 M Q in 
order to do so. An advantage of integrating the continuity equation 
backwards is that we will not have to make assumptions about the 
number density of BHs below 10 6 Mq. 



3.1 Deriving the X-ray luminosity function 

In order to learn whether the simple evolution scheme we de- 
scribed can explain observations and find possible values for the 
free parameters r), t\, t%, and Af const , we calculate the predicted 
hard XLF as a function of 2 and compare it to the luminosity- 
dependent density evo lution (LDDE) ha rd XLF derived from obser- 
vation in the study by Ueda et al. (2003). A numb er of subsequent 
studie s have also characte rized the XLF, includin g lLa Franca et alj 
d2005h.lAird et al.l fcOld) . and lUeda et all ( l201lh . Using the XLF 
lUeda et alj d2003h all ows o ne to more directly compare our results 
to iMerloni & He inz (2008), which uses the same XLF as a con- 
straint. 

The LDDE XLF of lUeda et all J2003I) is parametrized as fol- 
lows: 



d(t>(L 



x, 2) 



d log Lx 



A 



Lx 

u 



+ 



e(2,Lx) (9) 



with evolution term 



e(Z ' Lx) = ^[l + 2 c (Lx)r-(l + ^ *>Ze(L X ), (10) 



and redshift cut-off 



2 c (Lx) = 



L x > L a 



(11) 



The parameter values are: A = 5.04 x 10 Mpc ; L a = 



10 44t> erg s" 1 ; L t = 10* d a4 erg s" 1 ; 71 = 0.86; 72 = 2.23; 
pi = 4.23; p 2 - 1.5; z* = 1.9 and a = 0.335. 

In calculating the XLF from a given evolving BHMF, we will 
have to correct for the sampling bias of ob scured, Compton -thick 
AGN, which is not included in the XLF of lUeda et ail d2003h . The 
contribution of obscured sources is taken into account by divid- 
ing the BHMF by 1.6 before we use it to cal culate the XLF. This 
correc tion is similar t o the approach taken in Raimundo & Fabian! 
d2009h and lCao & LH d2008l) , and assumes the obscured source con- 
tribution to the XLF is independent of luminosity, although the pic- 
ture is likely to be more complicated and this simplifying assump- 
tion may be a re ason for discrepancie s we may find in the model 
and observations. [Treister et al.l (12010) find the ratio of obscured to 
unobscured quasars to be ~ 1 and increasing at redshifts 2 > 1.5. 

To calculate the XLF, we use the Eddington-ratio de- 
pendent bolometric correc tion Lx = Lboi/t(A) based on 
IVasudevan & Fabianl d2007» . We set: 



19.3 0.1 < A, 
«(A) = I 36.9 0.1 < A < 0.3, 
54.5 0.3 < A. 



(12) 



Assigning values for the free parameters 77, t\,t%, and M con st 
allows us to calculate (M{M, t)) and hence (m(M, t)). However, 
we will need to assume a distribution for m(M, t) in order to cal- 
culate the XLF. The spread in rh(M, t) is expected to be large for 
the following reasons. A significant fraction of BHs are known to 
be not accreting (or barely accreting) at a given time while oth- 
ers may be accreting highly. Also, the Eddington ratio A is a func- 
tion of m and a wide spread in Eddington ratios has been observed 
re.g. lKollmeieret all (12006) find a log-normal distribution of A for 
a fixed mass with a dispersion of ~ 0.3 dex in th e AGN and Galaxy 
Evolution Survey (AGES), and lRaimundo et al.l Q20 10) find a large 
number of low Eddington ratios 10~ 4 < A < 10 _1 in the Chan- 
dra Deep Fields]. It appears that the more sensitive an observa- 
tion is, the lower Eddington ratios may be found. The true spread 
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in m, or equiva l enfly, A, may even be larg e r tha n those found by 
iKollmeier et alj J200(J) or [kaimundo et al.l ( feoioh . since here we 
need to include all galaxies with supermassive BHs in the distri- 
bution, even ones that are not accreting. In this analysis, we will 
approximate BHs that are not accreting with BHs that have negli- 
gible accretion rh (their contribution to the total kinetic feedback 
will also be small). 

We assume a log-normal probability distribution for m(M, t). 
The distribution is centred around (rh(M,t)) — ^a 2 In 10 with 
standard deviation a (the centre is chosen so that the mean of 
the distribution is (m(M,t)), which is what we want). The stan- 
dard deviation a is the final free-parameter in our model. The log- 
normal distribution is a reasonable assumption because Edding- 
ton ratios appear to be roughly distributed in a lo g-normal fash- 
ion ( IKollmeier elafl 120061 : Raim undo et "all l20ld) so the m dis- 
tribution should be similar as well when we convert it using the 
relation between m and A. The shape of the left-tail of the distri- 
bution may not be well constrained by observations, but it has a 
small effect on the average of the distribution and objects in the left 
wing do not end up contributing to the 10 41 ' 5 -10 46 ' 5 erg s _1 XLF. 
Knowing the probability distribution D(rh(M,t)) of rh(M,t), 
one can obtain the Eddington ratio probability distribution using 
D(rh(M,t))drh = D(X(M,t))dX. To obtain the distribution of 
Lx, we separate the distribution D(X(M,t)) into three distribu- 
tions: 0.1 < A, 0.1 < A < 0.3, 0.3 < A (because of the dis- 
continuities in the bolometric correction factor k(A)), convert each 
using D(Lx(M,t))dLx = D(\(M,t))dX, and sum. Knowing 
the probability distribution of Lx for given M and t, as well as the 
number density of sources from the evolving BHMF at time t with 
masses between M and M + AM provides the information needed 
to construct the XLF. 

Calculating an XLF with the method described above for var- 
ious values of the free-parameters, we find that a does indeed have 
to be large (> 0.3 dex at least) in order to construct reasonable 
XLFs (the best-fitting value for all parameters is reported in §[4]l. 
Simply put, a < 0.3 dex does not produce objects with a wide 
enough distribution of Lx to cover at least 10 41 ' 5 to 10 46 ' 5 erg s _1 , 
which is the domain of the observationally derived XLF. The num- 
ber density of the BHs derived from an evolving BHMF is 1-2 
orders of magnitude higher than the number density of BHs in the 
observed XLF (for example, at z — 0, the number density of BHs 
is 2 x 10~ 2 Mpc~ 3 and the number density of X-ray sources with 
10 41 ' 5 < L x (erg s" 1 )" 1 < 10 46 5 is 3 x 10~ 4 Mpc -3 ). A large 
fraction of galaxies with BHs in fact do not contribute to the ob- 
served 10 41 ' 5 -10 46 ' 5 erg s _1 XLF, which is accounted for by hav- 
ing a large spread as dictated by a so that a considerable fraction of 
sources are accreting well below the critical accretion rate m cr it- 

Numerically calculating the XLF as a sum of Gaussians for 
given values of the free parameters rj, n, T2, Af cons t, and a, we 
find a good fit for the free parameters by minimizing the sum of the 



squared differences of log 
at log(L x (erg s" 



vcd( L X. z ) 



d log L x 



and log d< ^ dl g 



= 41.5,42.0,42.5, 



model is to be presented in §|4] 



log L x 

,46.5. The final 



3.2 Kinetic Feedback 

Once having found parameters for an evolutionary scenario that is 
able to reproduce observations (XLF and the local BHMF), the total 
amount of volume-averaged kinetic feedback may be calculated as 
a function of redshift, as well as a KLF, if we make an estimate for 
the fraction, /, of high accretion rate sources in HK mode (as op- 



posed to HR mode). This fraction / is the fraction of high-accreting 
sources that are radio-loud at a given cosmic time. We will present 
two calculations, one with a fiducial value of / = 0.1 across all 
redshifts and another with a refined estimate of / as a function of 
redshift needed to match the observed RLF for powerful, jetted ra- 
dio sources. 

The total kinetic energy per unit comoving volume at time t is 
calculated as: 



Lkin.tot (t) 



(A kin (M, t))L Edd (M) x dN i h f; t) di\ I (13) 



dM 



where 



(Akin (M,i))= / Akin (m)D(m(M,t)) dm. 



(14) 



The total kinetic energy per unit comoving volume may be broken 
down into the HK, HR, and LK contributions as well. 

The KLF is calculated similarly to the way it is done for the 
XLF. The distribution of rh(M, t) can be converted to a distribu- 
tion for Akm(Af , t) for both expressions for A k i n (Af, t) (due to the 
two high-accretion modes). The KLF may be calculated for both 
distributions and a weighted average may be taken according to the 
assumed value for / to obtain the total KLF. 



3.3 FR II sources and feedback 

From the calculated KLF, we can estimate the kinetic feedback due 
to powerful classical double-lobed FR II ( Fanaroff & Rilevl l 1974b 
sources and the amount of energy that goes into doing PdV work 
in expanding the lobes. Each lobe is powered by a jet with power 
Qjct = Lkin/2. We will consider sources with Q jct > 5 x 10 37 W 
to exhibit FR II j ets, which is the lower end of the dist ribution of jet 
power found in iBlundell. Rawli ngs & Willottl <fl999h for analytic 
models of FR II radio lobes to match observations. We can integrate 
the KLF (multiplied by Lkm) above luminosities Lkm = 10 38 W 
to find the total power going into FR II jets as a function of redshift. 
We estimate the fraction of the jet power that goes into PdV work 
to expand the lobes from the analyt i c mod el developed for FR II 
sources in lMocz. Fabian &Blundell (201 l|). The fraction is found 
by calculating J Q ' J pi dVi/(Qj e ttj) for a typical jet lifetime of tj = 
5 x 10 s yr and environmental para meters /3 = 1.5, ap = 10 kpc 
and 



1.67 x 10" 23 kg m" 3 ( IBlundell, Rawlings & Wiiiou] 
1999) and injection parameters 7 m j n = 1, 7max = 10 6 and p = 
2.14 |m ocz. F abian & Blundell 201 1). We find that the fraction, 
F, of power that goes into the expansion of the lobes for this set 
of parameters, as a function of jet power Qj at in units of W, can 
approximately be parametrized as: 



F(Qjet) = 2.4- 10 7 ( 



(15) 



For a typical jet power of Qj e t = 10 W, the fraction is ap- 
proximatel y 1/6. This value is reasonable compared to the cal- 
culation in Kaiser, De nnett-Thorpe & Alexanderl d 19971) (see their 
equation (17)) based on a simpler model of self-similar evolution 
of lobes. In our non-self-similar analytic model for evolution of 
FR II lobes, this fraction does depend somewhat on jet power and 
time along the evolution, and this dependence is being investigated 
in another paper (Mocz et al., in preparation). The remainder of 
the power that is supplied by the jet (other than expansion work 
by the lobes) goes into internal energy stored in the lobes, inter- 
nal energy in the hotspot (which is negligible), expansion work 
done by the hotspot, and energy lost due to synchrotron radiation 
in the radio and up-scattering of cosmic microwave background 
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(CMB) photons to X-ray energies via inverse-Comp ton (IC) scat- 
tering. We show in lMocz. Fabian & Blundelll ( l201ll) . that after the 
jets switch off, the radio emission due to synchrotron radiation 
plummets rapidly but the lobe may still be overpressured and will 
continue to expand, although not as rapidly as if the jet were still 
on. This means that there will still be PdV work done after the 
lobes are n ot detectable in the radio. Using the analytic model de- 
veloped in Mocz, Fabian & BlunddJ l l201ll) . we find that if a jet is 
turned on at z = 2 for tj = 5 x 10 8 yr, the ratio of PdV work 
the lobes do after the jets switch off until the present epoch z — 
is a considerable fraction, ~ 0.6, of the work done during the time 
the jet was on. The fraction of the total jet energy Qjetij lost due 
to synchrotron radiation is found to be negligible, ~ 0.002. The 
energy lost to inverse-Compton scattering of the CMB (ICCMB) 
may be more significant though still small, for example ~ 0.03 at 
2 = 2, but it is important to note that this fraction depends on red- 
shift due to the (1 + z) 4 dependence in the CMB energy density as 
a function of redshift. 

3.4 RLF for FR II sources 

In addition, with a KLF and a model for FR II radio sources, as 
well as an assumption about the jet duty cycles, we can derive an 
RLF from our popu lation and compare with the observed RLF of 
IWillott et al.l bOOll) . Comparison of the derived and observed RLF 
then tells us how to refine the jet duty cycle assumption / as a 
function of redshift. 

We assume that for a given redshift, the radio emission from 
an FR II source of jet power Q je t is randomly chosen between 
t = and 10tj in the evolution of the source (the radio emis- 
sion plummets very shortly after t = tj). The time t = 10tj 
corresponds to a time when the pressure in the lobes become 
comparable to the surrounding intergalactic medium (IGM; the 
lobes do remain overpressured for some time after jets switch off; 
iMocz. Fabian & BlundelilPioT l|). This is synonymous with saying 
that / = 0.1 of sources with potentially powerful jets have their jets 
turned on at a given time. The radio emissio n is calculated using the 
model in lMocz. Fabian & Blundelil 1201 lh . The RLF can then be 
derived from the KLF which tells us the number of sources with jet 
powers between Qj ct and Qj ct + AQj ct above Qj ct > 5 x 10 37 W 
(the minimum FR II jet power). 

The calculated KLF can be compared with the observed one 
and used to constrain / as a function of z. We can then eliminate 
our initial / = 0.1 assumption and derive more accurate an esti- 
mate for FR II feedback. 



4 RESULTS 

We carry out the steps mentioned in §[3] to evolve the BH popula- 
tion backwards in time and find that we can reasonably reproduce 
the observed XLF with our simple model. We find best-fitting pa- 
rameters r) = 0.08 ± 0.02, n = I0 9 ' 9±al yr, r 2 = 10 9 ' 9±01 yr, 
M const = 10 7 0±1 '° M , and a = 0.9 ± 0.1 dex, with 1-cr errors 
reported. The evolution of the BHMF is shown in Figure|2] A non- 
zero value for M CO nst is favoured, which leads to anti hierarchical 
growth, but the two time-scales n and T2 converge. The derived 
XLF from the population of BHs is presented in Figure [3] which 
we see agrees well with the observed XLF considering the simpli- 
fying assumptions we have made. The largest discrepancies (factor 
of 3 — 10) between the predicted and observed XLF occur at the 
'knee' feature of the observed XLF which is not possible to produce 



with the simple assumption of a smooth log-normal distribution of 
m and the low-luminosity end of the XLF, but this is generally 
within the 1-a Poisson errors in the observational determination 
of the XLF. A significant spread in the log-normal distribution of 
m, namely a = 0.9 ± 0.1 dex, is required to match observations 
well. The average values of m as a function of BH mass at various 
redshifts are plotted in Figure [4] The average values of m are well 
below 1 for these redshifts. 

Figure [5] shows the derived total kinetic feedback (and com- 
ponents due to different modes of accretion, and FR II sources) per 
unit comoving volume as a function of redshift assuming / = 0.1 
[meaning 1/10 of high-accretion sources have jets (i.e., are in HK 
mode) at a given cosmic time]. Figure [6] shows the accompanying 
derived KLF. These two figures are refined by then estimating an 
evolving duty cycle / = f(z) from deriving the RLF for FR II 
sources and comparing it to the observed RLF (Figure [7}. The evo- 
lution of the total kinetic feedback with this better informed duty 
cycle (plotted in Figure [8} is presented in Figure [9] and the cor- 
responding KLF is shown in Figure [TO) These are the three main 
figures of our paper. 

Our calculation of the RLF for FR II sources matches the ob- 
served RLF at z = 1 under the simple assumption that a fraction 
/ = 0.1 of powerful radio sources are switched on at a given time, 
corresponding to the picture that all sources have episodic jet ac- 
tivity with jets switched on for 5 x 10 8 yrs and enough time be- 
tween events so that the lobes stop being overpressured. With the 
/ = 0.1 assumption we overestimate the RLF at the present epoch, 
which means that a large fraction of once-powerful radio sources 
are switched off today, perhaps permanently if the fuel supply has 
depleted and nothing triggers more gas infall. At z > 1, we under- 
estimate the RLF, meaning that many radio sources have had to be 
switching on at redshifts 1.5 < z < 3 (the quasar era; this corre- 
sponds to a length of time of 4ij) and that a large fraction (> 0.1, 
close to 1 in fact) of high-accretion sources that could host power- 
ful jets were in fact active. Figure[8]shows, as a function of redshift, 
the fraction of high-accretion sources that need to be turned on for 
the derived and observed RLFs to match. The fraction varies from 
>0. 8 at 2 = 3 to 0.02 at 2 = 0. A revised estimate of the total ki- 
netic feedback is also shown in Figure [9] In such case, the amount 
of kinetic feedback rises to approximately match the radiative feed- 
back. The total kinetic feedback is then dominated by FR II sources 
at redshifts above z ~ 1.5. 



5 DISCUSSION 

We are able to reconstruct the observed XLF with our simple model 
for the growth of BHs. We integrate the local BHMF backwards 
according to a simple evolutionary scheme to obtain the evolving 
BHMF and average accretion rate as a function of BH mass and 
time. In our model, anti hierarchical BH growth is present due to a 
same quantity of mass, M CO nst, being accreted by all BHs irrespec- 
tive of BH mass. The term M cons t may be interpreted as a uniform 
source of matter available to be accreted by all BHs. Possibly such 
matter may be due to expelled material from an earlier, more vi- 
olent, merger-dominated accretion history of BHs before redshift 

Zl. 

The evolutionary scenario we described in this paper is simple 
with a limited number of free-parameters, which is not expected 
to fully describe in detail the evolutionary histories and feedback 
of individual supermassive BHs. However, the free-parameters all 
have distinct effects in the model, which avoids degeneracies and 
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Figure 2. The evolving BHMF that best reproduces the observed XLF. Anti 
hierarchical growth is present in the evolution. The alternating solid/dashed 
lines show the evolving BHMF from z = to z = 3 in increments of 
Az = 0.5. 
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Figure 3. Comparison of observed (dashed) and derived (solid) XLFs at 
z = (red), 1 (orange), 2 (green) and 3 (blue). The two are in relatively 
good agreement. Our simple assumption of a Gaussian spread in fa does not 
produce quite as sharp a 'knee' as in the observed XLF. Our simple model 
slightly over-predicts the XLF at z = 3. 

provides intuition on how the BHMF and XLF evolves. The free pa- 
rameters that describe the accretion history of supermassive BHs in 
our model are Mi, n, T2, and rj. The evolution of BHs in the high- 
mass end (M > Mi) is only primarily determined by 77 which ad- 
justs how much matter is accreting on to the BHs and n which de- 
scribes how this amount changes with time. BHs in the low-mass- 



10 

Log 10 (M BH / M ) 

Figure 4. The average value of m in our model at 2 = (red), 1 (orange), 2 
(green) and 3 (blue) as a function of BH mass, M. High and low mass BHs 
evolve differently (note the change in slope, due to the inclusion of M const 
in our evolution scenario) which allows for anti hierarchical growth. The 
parameters t\ and ti (which are found to be roughly the same) determine 
how quickly the low and high ends of the plot of fa changes with time. 

end (M < Mi) are allowed to evolve differently (different rates and 
change in rates from the high-mass counterparts) with the inclusion 
of the other two parameters Mi and T2 to allow for the possibility 
of anti hierarchical growth. Finally, we chose a log-normal distribu- 
tion of fn for a given mass characterized by a single free-parameter, 
the standard deviation a, to minimize degeneracies when match- 
ing the predicted and observed shape of the XLF. Supposing we 
have two different evolving BHMFs, it becomes impossible with 
a single parameter a to produce two XLFs with the same slopes 
and the same integral constraints. A more sophisticated model may 
improve the match between the observed and predicted XLF, but 
would be prone to degeneracies in the best-fitting free-parameters. 

Knowing the true distribution of fn for a given BH mass, 
which in general may change its shape and spread with redshift, 
is the most difficult to directly determine from the constraints we 
impose. We chose to investigate a simple log-normal distribution 
to understand what spread in the distribution the observations sug- 
gest. However, a refinement of the distribution would improve 
the predicted shape of the XLF function (although introduction 
of more than a single free parameter would introduce degenera- 
cies in the fitting of the XLF shape and hence larger uncertain- 
ties). A modified distribution for m, such as a log-normal with 
a power-law tail for low fn, a s suggested from observations by 
iKauffmann & Heckmarj (2009), could better reproduce the 'knee' 
feature in the XLF and the XLF at the lower luminosities (generally 
corresponding to lower-mass and smaller fn BHs). Such a modified 
distribution would mostly affect the number density and distribu- 
tion of low-rh sources, hence not affecting much our estimate of jet 
duty cycle with redshift, which is determined by the evolution of 
high-ra sources. 

The large spread we find in m, namely a standard deviation of 
a = 0.9 ± 0. 1 dex, helps account for BHs that are not or negligibly 
accreting. This spread corresponds to a ~ 0.9 dex standard devi- 
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Figure 5. The energies of various components of BH feedback per comov- 
ing volume as a function of redshift. Contribution due to only unobscured 
sources is shown; to estimate the effect of the obscured sources, multiply by 
1.6. The total radiative feedback (total bolometric luminosity of sources) is 
indicated by the thick red line. The light blue shaded region shows the pos- 
sible total kinetic feedback if one varies the fraction / of sources accreting 
above rh cr ; t in HK mode from / = to 1. The thick dashed blue line 
shows the kinetic output for / = 0.1. The orange, green, and purple lines 
show the breakdown contribution of LK, HK and HR sources, respectively, 
to the total kinetic feedback. The kinetic feedback due to FR II sources 
(^kin = 2Qj et > 10 38 erg s _1 ) is shown in dashed blue and closely 
follows the HK line. The fraction of the FR II kinetic power that goes into 
PdV work to expand the lobes is shown in dotted blue. 



ation spread in Eddington ratios. A significant fraction of sources 
thus have A < 10 , corresponding to sources with little/negligible 
accretion activity and sources that do not contribute to the XLF. As 
the average accretion rate, (rh(M,t)), decreases with time in our 
model, a higher fraction of sources have little activity (A < 10 -3 ) 
as we reach z = 0. 

From our evolution scenario we can quantify the amounts of 
different modes of feedback and learn which one dominates. The 
amount of kinetic feedback may exceed radiative feedback in two 
cases. First, if the fraction of HK sources is high, then the kinetic 
feedback will always exceed radiative feedback at all redshifts (if 
we assume the BH accretion physical model that is presented Fig- 
ure [TJ. Secondly, once most of the sources are accreting in LK 
mode, kinetic feedback will exceed radiative feedback, which oc- 
curs in our calculations for z < 0.4. In our calculations, with the 
refined assumption of the duty cycle / changing throughout cos- 
mic time, the total kinetic and radiative feedback closely match 
each other, which is somewhat of a coincidence caused by both the 
average accretion rates and duty cycle / decreasing with smaller 
redshift. 

The total kinetic feedback may be decomposed into HK, HR, 
and LK mode feedback. The decomposition is important to under- 
stand because different fractions of the energy from the three feed- 
back modes may go into actually heating the interstellar medium 
(ISM) and IGM. We see from Figure [9] that, with our estimate for 
evolving jet duty cycle, much of the total kinetic output energy at 
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Figure 6. The derived KLF from our simple evolutionary scenario, assum- 
ing a fraction / = 0.1 of highly accreting sources are HK and a fraction of 

1 — / = 0.9 highly accreting sources are HR. The KLF is shown at z = 
(red), 1 (orange), 2 (green) and 3 (blue). The dashed lines correspond to 
the HK plus / = 0.1 of LK sources and the dotted lines correspond to 
the HR plus 1 — / = 0.9 of LK sources. Sources with kinetic luminosi- 
ties .Lkin > 10 45 erg s" 1 (corresponding to the minimum jet energy we 
consider to create FR II type lobes) are mostly HK sources. 

2 < 1.5 comes from the LK mode, and its contribution to the total 
power may increase up to z ~ 0.8 even though the total kinetic 
feedback decreases with time. Our calculations highlight the im- 
portance that sources in LK mode may have in feedback. Since the 
total feedback in LK mode is higher than HK and HR modes at 
z < 1.5, even if a smaller fraction of the LK mode kinetic en- 
ergy goes into heating the gas in the host galaxy it could still be 
more important. Sour ces accreting in L K mode would correspond 
to FR I radio sources. |Pe Yound f2010) argued, using comparisons 
with observed data, that feedback in FR I sources may occur due 
to AGN outflows that are strongly decelerated and become fully 
turbulent sonic or subsonic flows due to their interaction with the 
surrounding medium. 

Our calculations also show, however, that at earlier redshifts 
2 > 1.5, when the jet duty cycle is large (> 0.5), the total volume- 
averaged energy feedback from the HK mode (classical double- 
lobed FR II type feedback) is dominant. These sources expand 10- 
1000 kpc lobes and do PdV work on the surrounding IGM. The 
amount of PdV work during the time the jet is on is found to be 
roughly 1 /6 of the kinetic energy that goes into the jets. It is un- 
clear how much of this energy goes into heating the central galaxy, 
which can be presumably small since the hotspots and lobes are far 
from the galaxy for the majority of the time the jet is switched on, 
and rather these giant lobes give considerable feedback to the IGM. 
Galaxies that show FR II processes must have had feedback heat- 
ing their ISM as well, and perhaps much of this feedback comes 
either from radiative feedback or the kinetic feedback during the 
time spent in HR mode. If the HR mode is the predominant form of 
feedback in the ISM for FR II type sources, then a much larger frac- 
tion of the output energy is needed to heat the ISM than would be 
needed if HK mode were the culprit, as the total volume- averaged 
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Figure 7. The derived RLF for FR lis from the KLF. The thick dashed 
lines correspond to the derived RLFs at z = (red), 2 = 1 (orange), 
2 = 2 (green) and 2 = 3 (blue) assuming that / = 0.1 of such sources 
are switched on at all redshifts. The thin li nes are the observed RLFs for 
the population of the brightest FR lis from Willott et al There is 

discrepancy between the two sets of RLFs, due to the / = 0.1 assumption. 
Underestimation means that > 0.1 potential FR II sources are switched on 
at that redshift. Overestimation means that < 0.1 sources are on. The RLF 
is corrected by varying / with redshift and the corrected RLF is shown in 
thick solid lines. 




Figure 8. Approximate fraction of powerful radio sources switched on (jet 
duty cycle) as a function of redshift so that our derived RLF matches the 
observed one. Dashed line shows a fiducial 0.1 value for comparison. This 
amounts to correcting the solid kinetic line (/ = 0.1) in Fig.[5]to the ob- 
served (dot-dashed) kinetic. 
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Figure 9. Same as Figure \5\ except with the refined assumption that / can 
vary with redshift, which is estimated from the observed RLF. It is interest- 
ing to note that the total kinetic output follows the radiative energy output 
closely. Also, with this refined estimate for / we see that HK mode feed- 
back becomes much more dominant at 2 > 1.5. 
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Figure 10. Same as Figure |6]except with the refined assumption that / can 
vary with redshift, which is estimated from the observed RLF. 



energy from HR feedback is an order of magnitude below HK feed- 
back at z > 1.5. 

With our approach to determine the total kinetic feedback 
(which uses the RLF along with a model for FR II evolution 
to estimate the jet-duty cycle as a function of redshift) we find 
a significantly higher estimate than previous studies by almost 
an order of magnitude at z > 1.5. We predict that the to- 
tal volume averaged kinetic feedback closely traces the radia- 



© 2013 RAS, MNRAS 000,[T]03] 



10 P. Mocz et. al. 



tive fe edback rather than bein g about an order of magnitude be- 
low it. ICattaneo & Bestl d2009h also calculate radiative and kinetic 
feedback as a function of redshift, u sing the bolometric LF of 
Hopki ns, Richards & Hernquist ( 2007) to find radiative feedback 
and the RLF of iDunlop & Peacock! q 199CT) for extragalactic radio 
sources to estimate mechanical power in kinetic feedback. Their 
calculated ratio for radiative to kinetic feedbac k agrees with the 
earlier work of iKdrding. Jester & Fenden d2008l) . who also use an 
RLF directly for calculating the mechanical feedback (unlike our 
work). We find a higher result for total kinetic feedback because 
some LK sources may be missing from the RLF, as well as a frac- 
tion of br ight sources due to the 'you th-redshift' degeneracy de- 
scribed in iBlundell & Rawlingsl ( fl999h . which says that older radio 
sources, as they expand, decrease their synchrotron radiation in the 
radio (due to lowering of the magnetic field strengths during expan- 
sion) and fall below the flux limit of the surveys used to construct 
the RLF and thus go undetected. A model for the time evolution of 
the radio emission of radio sources is needed to accurately have a 
handle on the 'youth-redshift' degeneracy, which is part of our ap- 
proach. Our calculation of radiative feedback, u nlike total kinetic 
feedba ck, agrees very well with the calculation in lCattaneo & Bestl 
(2009). The amount of PdV work by FR II sources is comparable 
to the amount of cavity work calculated in lCattaneo & Bestl (2009) 
as a function of redshift. But we find that the total kinetic feedback 
is significantly higher than the fraction that goes into expanding 
lobes, and is comparable to the total radiative feedback. 

Kinetic radio-mode feedback is previously overlooked and po- 
tentially very important for the grow th and feedbac k in galaxies on 
a cosmological scale. The review by |Fabiart ( l2012T) identifies FR II 
feedback (powerful, jetted radio outbursts) as an ill-understood 
feedback mode which can deposit the kind of energy that can dis- 
rupt even a group of galaxies. In comparison to FR I kinetic feed- 
back and radiative (quasar) feedback, FR II feedback has been little 
discussed in the literature. Here we have provided a first estimate of 
the amplitude of FR II feedback. From Figure|9]we see that its am- 
plitude may be very significant above z ~ 1.5, rivalling radiative 
output, much of which passes beyond the host galaxy. Accounting 
for this level of FR II type feedback in cosmological simulations 
may refine our picture of galaxy evolution. 
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